<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Minimize transition bandwidth of a linear phase lowpass FIR filter</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/filter_design/html/fir_lin_phase_lowpass_min_trans.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Minimize transition bandwidth of a linear phase lowpass FIR filter</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs a linear phase FIR lowpass filter such that it:</span>
<span class="comment">% - minimizes the transition band width (i.e. minimize w_stop)</span>
<span class="comment">% - has a constraint on the maximum passband ripple</span>
<span class="comment">% - has a constraint on the maximum stopband attenuation</span>
<span class="comment">%</span>
<span class="comment">% This is a quasiconvex problem and is solved using a bisection.</span>
<span class="comment">%</span>
<span class="comment">%   minimize   w_stop</span>
<span class="comment">%       s.t.   1/delta &lt;= H(w) &lt;= delta     for w in the passband</span>
<span class="comment">%              |H(w)| &lt;= atten_level        for w in the stopband</span>
<span class="comment">%</span>
<span class="comment">% where H is the frequency response function and variable is</span>
<span class="comment">% the filter impulse response h (and its order/length).</span>
<span class="comment">% Data is delta (max passband ripple) and atten_level (max stopband</span>
<span class="comment">% attenuation level).</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% user's filter specifications</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% starting point for the stopband (needs to be feasible)</span>
wstop = 0.24*pi;        <span class="comment">% stopband start freq (in radians)</span>
TOL = 1e-3;             <span class="comment">% precision to which we should run bisection</span>

n = 10;                 <span class="comment">% filter order (2n+1 is the full order)</span>
wpass = 0.12*pi;        <span class="comment">% passband cutoff freq (in radians)</span>
delta = 1;              <span class="comment">% max (+/-) passband ripple in dB</span>
atten_level = -30;      <span class="comment">% stopband attenuation level in dB</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% create optimization parameters</span>
<span class="comment">%********************************************************************</span>
m = 30*n; <span class="comment">% freq samples (rule-of-thumb)</span>
w = linspace(0,pi,m);

<span class="comment">%*********************************************************************</span>
<span class="comment">% use bisection algorithm to solve the problem</span>
<span class="comment">%*********************************************************************</span>

wstop_bot = wpass;
wstop_top  = wstop;

<span class="keyword">while</span>( wstop_top - wstop_bot &gt; TOL)
  <span class="comment">% try to find a feasible design for given specs</span>
  wstop_cur = (wstop_top + wstop_bot)/2;

  <span class="comment">% create optimization matrices (matrix of cosines)</span>
  A = [ones(m,1) 2*cos(kron(w',[1:n]))];

  <span class="comment">% passband 0 &lt;= w &lt;= w_pass</span>
  ind = find((0 &lt;= w) &amp; (w &lt;= wpass));    <span class="comment">% passband</span>
  Ap  = A(ind,:);

  <span class="comment">% transition band is not constrained (w_pass &lt;= w &lt;= w_stop)</span>

  <span class="comment">% stopband (w_stop &lt;= w) (this is the changing constraint)</span>
  ind = find((wstop_cur &lt;= w) &amp; (w &lt;= pi));   <span class="comment">% stopband</span>
  As  = A(ind,:);

  <span class="comment">% formulate and solve the feasibility linear-phase lp filter design</span>
  cvx_begin <span class="string">quiet</span>
    variable <span class="string">h_cur(n+1,1)</span>;
    <span class="comment">% feasibility problem</span>
    <span class="comment">% passband bounds</span>
    Ap*h_cur &lt;= 10^(delta/20);
    Ap*h_cur &gt;= 10^(-delta/20);
    <span class="comment">% stopband bounds</span>
    abs( As*h_cur ) &lt;= 10^(atten_level/20);
  cvx_end

  <span class="comment">% bisection</span>
  <span class="keyword">if</span> strfind(cvx_status,<span class="string">'Solved'</span>) <span class="comment">% feasible</span>
    fprintf(1,<span class="string">'Problem is feasible for stopband freq = %3.4f rads\n'</span>,wstop_cur);
    wstop_top = wstop_cur;
    <span class="comment">% construct the full impulse response</span>
    h = [flipud(h_cur(2:end)); h_cur];
  <span class="keyword">else</span> <span class="comment">% not feasible</span>
    fprintf(1,<span class="string">'Problem is not feasible for stopband freq = %3.4f rads\n'</span>,wstop_cur);
    wstop_bot = wstop_cur;
  <span class="keyword">end</span>
<span class="keyword">end</span>

wstop = wstop_top;
fprintf(1,[<span class="string">'\nOptimum stopband frequency for given specs is %3.4f*pi rads\n'</span> <span class="keyword">...</span>
           <span class="string">'and the minimum transition width is %3.4f*pi radians.\n'</span>],<span class="keyword">...</span>
            wstop/pi, (wstop-wpass)/pi);


<span class="comment">%********************************************************************</span>
<span class="comment">% plots</span>
<span class="comment">%********************************************************************</span>
figure(1)
<span class="comment">% FIR impulse response</span>
plot([0:2*n],h',<span class="string">'o'</span>,[0:2*n],h',<span class="string">'b:'</span>)
xlabel(<span class="string">'t'</span>), ylabel(<span class="string">'h(t)'</span>)

figure(2)
<span class="comment">% frequency response</span>
H = exp(-j*kron(w',[0:2*n]))*h;
<span class="comment">% magnitude</span>
subplot(2,1,1)
plot(w,20*log10(abs(H)),<span class="keyword">...</span>
     [wstop pi],[atten_level atten_level],<span class="string">'r--'</span>,<span class="keyword">...</span>
     [0 wpass],[delta delta],<span class="string">'r--'</span>,<span class="keyword">...</span>
     [0 wpass],[-delta -delta],<span class="string">'r--'</span>);
axis([0,pi,-40,10])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'mag H(w) in dB'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
Problem is not feasible for stopband freq = 0.5655 rads
Problem is not feasible for stopband freq = 0.6597 rads
Problem is feasible for stopband freq = 0.7069 rads
Problem is not feasible for stopband freq = 0.6833 rads
Problem is feasible for stopband freq = 0.6951 rads
Problem is not feasible for stopband freq = 0.6892 rads
Problem is not feasible for stopband freq = 0.6921 rads
Problem is feasible for stopband freq = 0.6936 rads
Problem is not feasible for stopband freq = 0.6929 rads

Optimum stopband frequency for given specs is 0.2208*pi rads
and the minimum transition width is 0.1008*pi radians.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_lin_phase_lowpass_min_trans__01.png" alt=""> <img src="fir_lin_phase_lowpass_min_trans__02.png" alt=""> 
</div>
</div>
</body>
</html>